c     C  O N G P 2 D 3   last revised Feb 6, 1992
c     This program is made available courtesy of Rose Byrne
c     Feb 10, 1992.  Rose is a doctoral candidate in Civil Engineering
c     at UVa.
c
c     Congp2d3.for is a revision of congp2dr.for.  Congp2d3.for
c     extends the algorithm described in Roache, "Computational
c     Fluid Mechanics", to an irregular, simply connected,
c     possibly non-convex region.  If the region is non-convex, the
c     indentation should be placed in such a way that there is no
c     node j which is neither physically adjacent to node j-1
c     nor along the side whose nodes are specified in nbot.
c     The algorithm in Roache originally appeared in Sundberg,
c     "Two Computerized Methods For Plotting Functions of Two
c     Independent Variables". This algorithm involves searching each
c     element for sides which contain points on a given contour.
c     The points within the element are then connected, and the pen
c     is lifted.  Even complicated results can be contoured
c     WITHOUT MANUALLY REORDERING THE DATA POINTS.  In the rare cases
c     when a contour intersects a cell in four places, the program
c     prepares a data file which connects those four dots in all 6
c     possible ways.  The user should plot the file first and determine
c     from the context of the plot how those four points should be
c     connected.  Inserting a blank line in the data file will cause
c     gnuplot to pick up the pen.  Congp2d3.for finds the points and
c     writes a data file for gnuplot with blank lines in each place
c     where the pen should be lifted.  Congp2d3.for also lifts the pen
c     in each place where the region is non-convex.  In each place
c     where the region is non-convex, the contours must stop at the
c     boundary of the region and pick up on the other side.
c     Sundberg's algorithm as extended here treats contours on
c     non-convex regions correctly without the user's intervention.
c     However, the section of the code which identifies points on the
c     boundary cannot distinguish between a rectangular indentation
c     in the region and a trapezoidal indentation.  The points which
c     are on the boundary of the region but are also in the same
c     column as another point on the boundary will not be found and must
c     be entered manually.  For most regions, this will be only a couple
c     of points whose coordinates will be obvious.  Process all data
c     files for a given region and then fix the data file for the
c     boundary.
c
c
c     Congp2d3.for assumes that the points to be contoured are
c     numbered as they usually are for finite element methods:
c     from the bottom to the top of each column, with each column
c     beginning where the last one leaves off.  Congp2d3.for keeps
c     track of which nodes are in each element using an array of
c     the points at the bottom of the region.  The following
c     is an example of how the points are numbered:
c
c           5-----10               23----28
c           |     |                 |     |    33
c           |     |                 |     |     |
c           4-----9----14----18----22----27----32----37
c           |     |     |     |     |     |     |     |    41
c           |     |     |     |     |     |     |     |     |
c           3-----8----13----17----21----26----31----36----40
c           |     |     |     |     |     |     |     |     |   44
c           |     |     |     |     |     |     |     |     |    |
c           2-----7----12----16----20----25----30----35----39---43
c           |     |     |     |     |     |     |     |     |    | 46
c           |     |     |     |     |     |     |     |     |    |  |
c           1-----6----11----15----19----24----29----34----38---42-45
c
c           the nodes on the bottom are 1,6,11,15,19,24,29,34,38,42,45
c
c     This program processes a certain kind of data file which cannot
c     be contoured by gnuplot3.0 (even after re-ordering the points) so
c     that gnuplot can construct a contour plot by drawing the contour
c     lines and domain boundary on a 2-D grid using the 2-D version
c     of gnuplot.  This program processes a data file such as that
c     produced by a finite element or finite difference program which
c     contains values of the dependent variable at each point of a 2-D
c     possibly irregular grid on a possibly irregular domain.  The
c     elements are not necessarily rectangular.  If triangular elements
c     are used in the finite element model which generates the data,
c     the side of a triangular element which is on the boundary will not
c     be searched for contour crossings, but the elements on the
c     interior of the region will be combined to make quadrilateral
c     elements.  Non-rectangular quadrilaterals work fine.
c
c     This program reads in the three arrays containing the three
c     coordinates (dependent variable and two independent variables)
c     of points representing the data which the user would
c     like to have contoured.  These points may correspond to the points
c     of any finite element or finite difference grid.  The arrays are
c     assumed to be ordered in the order that the grid points are
c     numbered. It is very helpful if the order in which the grid points
c     are numbered is nowhere near parallel to the expected direction of
c     the contour lines.  The program uses linear interpolation to
c     construct level curves at equally spaced intervals:  the
c     curves which will form the contours.  The coordinates of these
c     curves are written to a file in the form which gnuplot requires.
c     The curves will not contain the same number of points, so a
c     contour plot cannot be constructed using the contour option of
c     the 3-D version of gnuplot.  A contour plot can, however, be
c     constructed using the 2-D version of gnuplot by drawing the
c     independent variable coordinates of each contour line on the
c     same 2-D grid.
c
      real zero
      integer npts,nel,nlvl,lbot,mxlen
c
c     define parameters so dimension statements don't have to be
c     changed each time the data is changed, yet arrays are no
c     larger than they need to be.  The parameter statements must be
c     changed each time the 2-D grid to be contoured or the number of
c     contours desired is changed.
c     npts = the number of points in the 2-D grid.
c     nel = the number of elements in the 2-D grid.
c     nlvl = the number of equally spaced level curves (contours)
c            to be drawn.
c     lbot = the number of points along the bottom (or lhs) of the grid.
c     mxlen = the maximum number of points allowed in a contour.
      parameter(npts=595)
      parameter(nel=528)
      parameter(nlvl=15)
      parameter(lbot=52)
      parameter(mxlen=400)
c
c     The program uses the following arrays:
c     x = the array of the independent variable in one direction.
c     z = the array of the independent variable in the other direction.
c     h = the array of the dependent variable.
c     xcont = the array of x coordinates of points on the contours.
c            The contents, and the contents of zcont and the
c            appropriate value of h for each contour will be
c            written to the output file for processing by gnuplot.
c     zcont = the array of z coordinates of points on the contours.
c     len = the number of contour points in each contour line, not
c           counting repetitions for Sundberg's connect the dots
c           algorithm but counting each dot in each element in which
c           it is found.
c     value = the value of the dependent variable along a given contour.
c     nbot = the array of points along the bottom of the grid (if points
c            are numbered bottom to top).  If points are numbered left
c            to right, this array contains the points on the left-hand
c            edge. npts+1 is always added to the end of the array
c     ncell = the number of contour points found in each element
c     last = a running tally of the number of points in this contour
c            level including repetitions for Sundber's algorithm.
      real x(npts),z(npts),h(npts),xcont(nlvl,mxlen)
      real value(nlvl), zcont(nlvl,mxlen)
      integer ncell(nlvl,nel), last(nlvl), len(nlvl), nbot(0:lbot)
c
c     call subroutine readin to read from the screen the names of the
c     input and output files and variables not included in the files,
c     to open the necessary files, and to read the data from the input
c     files.
      call readin(zero,npts,nlvl,lbot,h,x,z,nbot,value,len)
c
c     call subroutine scan to scan all points in the array, looking for
c     places where a contour value should be.
      call scan(h,x,z,nbot,value,len,ncell,xcont,zcont,
     $     npts,nel,nlvl,lbot,mxlen,zero)
c
c     call subroutine wrcon to write the information about the contours
c     to a file in the correct format for gnuplot.
      call wrcon(ncell,xcont,zcont,last,nlvl,nel,mxlen)
c
c     call subroutine wrbdy to write the information about the
c     coordinates of the points on the boundary to another file
c     for use by gnuplot.
      call wrbdy(x,z,nbot,npts,lbot)
c
c     call subroutine wrlabl to write the labels for the contours
c     to a third file for use by gnuplot.
      call wrlabl(xcont,zcont,value,len,nlvl,mxlen)
c
c     To plot the data using gnuplot, first enter gnuplot by typing
c     gnuplot with the dafault drive set to the drive that gnuplot
c     is on.  Set any desired titles, axis labels, and output device
c     (with appropriate options).  See the gnuplot documentation for
c     details.  Then type the following commands:
c         load "labels"
c         plot "outcon", "outbnd"
c     Insert the appropriate file names for labels, outcon, and outbnd
c     above.  The quotation marks are part of the commands.
c     See the gnuplot documentation for how to set axis ranges and
c     many other features of the plot.
      stop
      end
c
      subroutine readin(zero,npts,nlvl,lbot,h,x,z,nbot,
     $     value,len)
      character*16 incont,outcon,rdnbot,outbnd,labels
      integer npts, nlvl, lbot, iflag
      real zero, to, from, by
      real h(npts), x(npts), z(npts), value(nlvl)
      integer nbot(0:lbot), len(nlvl)
c     prompt user for input and output file names
      write(*,*) 'Welcome to congp2d3, a preprocessor for gnuplot.'
      write(*,*) 'This program will take the output from your 2-D'
      write(*,*) 'finite element or finite difference program'
      write(*,*) '(or a similar data file) and put it in a form'
      write(*,*) 'in which it can be used by gnuplot.'
      write(*,*)
      write(*,*) 'Please type the name of the file with your results'
      read(*,'(a16)') incont
      write(*,*) 'Please type the name of the file with the array of'
      write(*,*) 'points along the bottom of the grid (for points'
      write(*,*) 'numbered bottom to top'
      read(*,'(a16)') rdnbot
      write(*,*) 'Thank you.  Now please type the name of the file to'
      write(*,*) 'store the contour data for use by gnuplot.'
      read(*,'(a16)') outcon
      write(*,*) 'Thank you.  Now please type the name of the file to'
      write(*,*) 'store the boundary data for use by gnuplot.'
      read(*,'(a16)') outbnd
      write(*,*) 'Thank you.  Now please type the name of the file to'
      write(*,*) 'store the label definitions for gnuplot'
      read(*,'(a16)') labels
      write(*,*) 'Thank you.  Now please type a positive number which'
      write(*,*) 'is close enough to zero (for use in identifying'
      write(*,*) 'points on the boundary which should be on one of'
      write(*,*) 'the contours).'
      read(*,*) zero
      write(*,*) 'Thank you.'
      write(*,*) 'The parameter statements indicate that ',nlvl
      write(*,*) 'contours will be drawn.  Please type the value of'
      write(*,*) 'the dependent variable along the contour where that'
      write(*,*) 'dependent variable is maximum.'
      read(*,*) to
      write(*,*) 'Thank you.  Please type the value of the dependent'
      write(*,*) 'variable along the contour where that dependent'
      write(*,*) 'variable is minimum.'
      read(*,*) from
      by = (to - from)/float(nlvl-1)
      write(*,*) 'Thank you.  I compute the contour interval as',by
      write(*,*) 'If something is wrong, please type 0 to stop the'
      write(*,*) 'program.  Then edit the parameter statements.'
      write(*,*) 'Otherwise, type 1 to continue with the program.'
      read(*,*) iflag
      if(iflag .eq. 0) stop
c
c     open files with variable file names to save wear and tear on user
      open(9,file=incont)
      open(10,file=outcon)
      open(11,file=rdnbot)
      open(12,file=outbnd)
      open(13,file=labels)
c
c     read in the arrays x,z,h,nbot
      read(9,*) (h(j),j=1,npts)
      read(9,*) (x(j),j=1,npts)
      read(9,*) (z(j),j=1,npts)
      read(11,*) (nbot(jcol),jcol=0,lbot)
c
c     find the values of the dependent variable along the contours
      write(*,*) 'values of the dependent variable along the contours'
      do 100 i=1,nlvl
           value(i) = from + (i-1)*by
           len(i) = 0
           write(*,*) value(i)
  100 continue
      return
      end
c
      subroutine scan (h,x,z,nbot,value,len,ncell,xcont,zcont,npts,
     $     nel,nlvl,lbot,mxlen,zero)
      integer npts, nel,nlvl,lbot,mxlen
      real h(npts), x(npts), z(npts), xcont(nlvl, mxlen),
     $     zcont(nlvl,mxlen), value(nlvl)
      integer nbot(0:lbot), len(nlvl), ncell(nlvl,nel)
c     for each contour, scan all points in the array, looking for
c     adjacent points where the dependent variable is bigger than
c     the contour value at one point and smaller than the contour
c     value at the other point.  Let the current point be a candidate
c     for the upper lefthand corner of some element.  If the lower
c     righthand corner and the upper righthand corner of that element
c     exist, check for contour values between the upper lefthand corner
c     and the lower lefthand corner, between the upper lefthand corner
c     and the upper righthand corner, between the lower lefthand corner
c     and the lower righthand corner, and between the upper righthand
c     corner and the lower righthand corner.  If the lower righthand
c     corner or the upper righthand corner does not exist, any needed
c     contour points have already been found. Keep track of how many
c     points are found for each element.  The appropriate element of
c     ncell is the number of contour points in this cell for this
c     contour level.  i counts the contour levels and icell counts the
c     points in the contour level.
      do 500 i=1,nlvl
           val = value(i)
           icell = 0
           do 400 jcol = 0,lbot-2
                nbottm = nbot(jcol)
                nxtbot = nbot(jcol+1)
                lcol = nxtbot - nbottm
                ntop = nxtbot - 1
                ntopn = nbot(jcol+2) - 1
                do 300 j = nbottm + 1, ntop
c                 do lower righthand corner and upper righthand corner
c                 of this element exist?  If not, go to the next
c                 element.
                  jp = j+lcol
                  if (jp .le. ntopn .and. jp-1 .le. ntopn) then
                     icell = icell + 1
                     ncell(i,icell) = 0
c                    check for a contour point between the upper
c                    lefthand corner and the lower lefthand corner
                     if ( (val .lt. h(j) .and. val .gt. h(j-1)) .or.
     $                    (val .gt. h(j) .and. val .lt. h(j-1)) ) then
c                         the contour passes through a point between
c                         points j and j-1.  Find it by linear
c                         interpolation.
                          frac = (h(j) - val) /(h(j) - h(j-1))
                          dx = x(j) - x(j-1)
                          dz = z(j) - z(j-1)
                          xval = x(j) - frac*dx
                          zval = z(j) - frac*dz
c                         xval,zval is the point to plot.
                          len(i) = len(i) + 1
                          xcont(i,len(i)) = xval
                          zcont(i,len(i)) = zval
                          ncell(i,icell) = ncell(i,icell) + 1
                     end if
                     jp = j + lcol
c                    check for a contour point between the upper
c                    lefthand corner and the upper righthand corner.
                     if ( (val .lt. h(j) .and. val .gt. h(jp)) .or.
     $                      (val .gt. h(j) .and. val .lt. h(jp)) ) then
c                           the contour passes through a point between
c                           points j and jp.  Find it by linear
c                           interpolation.
                            frac = (h(j) - val) /(h(j) - h(jp))
                            dx = x(j) - x(jp)
                            dz = z(j) - z(jp)
                            xval = x(j) - frac*dx
                            zval = z(j) - frac*dz
c                           xval,zval is the point to plot.
                            len(i) = len(i) + 1
                            xcont(i,len(i)) = xval
                            zcont(i,len(i)) = zval
                            ncell(i,icell) = ncell(i,icell) + 1
                     end if
c                    check for a contour point between the upper
c                    righthand corner and the lower righthand corner.
                     if ( (val .lt. h(jp) .and. val .gt. h(jp-1)) .or.
     $                    (val .gt. h(jp) .and. val .lt. h(jp-1)) ) then
c                         the contour passes through a point between
c                         points jp and jp-1.  Find it by linear
c                         interpolation.
                          frac = (h(jp) - val) /(h(jp) - h(jp-1))
                          dx = x(jp) - x(jp-1)
                          dz = z(jp) - z(jp-1)
                          xval = x(jp) - frac*dx
                          zval = z(jp) - frac*dz
c                         xval,zval is the point to plot.
                          len(i) = len(i) + 1
                          xcont(i,len(i)) = xval
                          zcont(i,len(i)) = zval
                          ncell(i,icell) = ncell(i,icell) + 1
                     end if
c                    check for a contour point between the lower
c                    lefthand corner and the lower righthand corner.
                     if ( (val .lt. h(j-1) .and. val .gt. h(jp-1)) .or.
     $                      (val .gt. h(j-1) .and. val .lt. h(jp-1)) )
     $                      then
c                           the contour passes through a point between
c                           points j-1 and jp-1.  Find it by linear
c                           interpolation.
                            frac = (h(j-1) - val) /(h(j-1) - h(jp-1))
                            dx = x(j-1) - x(jp-1)
                            dz = z(j-1) - z(jp-1)
                            xval = x(j-1) - frac*dx
                            zval = z(j-1) - frac*dz
c                           xval,zval is the point to plot.
                            len(i) = len(i) + 1
                            xcont(i,len(i)) = xval
                            zcont(i,len(i)) = zval
                            ncell(i,icell) = ncell(i,icell) + 1
                     end if
                     if (abs(val - h(j)) .le. zero) then
c                         check for a contour point at j
c                         (upper lefthand corner)
c                         check for val being equal to h(j)
                          len(i) = len(i) + 1
                          xcont(i,len(i)) = x(j)
                          zcont(i,len(i)) = z(j)
                          ncell(i,icell) = ncell(i,icell) + 1
                     end if
                     if (abs(val - h(j-1)) .le. zero) then
c                         check for a contour point at j-1
c                         (lower lefthand corner)
c                         check for val being equal to h(j-1)
                          len(i) = len(i) + 1
                          xcont(i,len(i)) = x(j-1)
                          zcont(i,len(i)) = z(j-1)
                          ncell(i,icell) = ncell(i,icell) + 1
                     end if
                     if (abs(val - h(jp)) .le. zero) then
c                         check for a contour point at jp
c                         (upper righthand corner)
c                         check for val being equal to h(jp)
                          len(i) = len(i) + 1
                          xcont(i,len(i)) = x(jp)
                          zcont(i,len(i)) = z(jp)
                          ncell(i,icell) = ncell(i,icell) + 1
                     end if
                     if (abs(val - h(jp-1)) .le. zero) then
c                         check for a contour point at jp-1
c                         (lower righthand corner)
c                         check for val being equal to h(jp-1)
                          len(i) = len(i) + 1
                          xcont(i,len(i)) = x(jp-1)
                          zcont(i,len(i)) = z(jp-1)
                          ncell(i,icell) = ncell(i,icell) + 1
                     end if
                  end if
  300           continue
  400      continue
  500 continue
      return
      end
c
      subroutine wrcon(ncell,xcont,zcont,last,nlvl,nel,mxlen)
      integer nlvl,nel,mxlen
      real xcont(nlvl,mxlen), zcont(nlvl,mxlen)
      integer ncell(nlvl,nel), last(nlvl)
c     write information needed by gnuplot to file in format appropriate
c     for gnuplot.  Each blank line causes gnuplot to lift the pen.
c     Write the contour points for each cell with the order and
c     repetition required for Sundberg's algorithm.  Then lift the
c     pen and go to the next cell.
      do 700 i=1,nlvl
           last(i) = 0
           do 650 icell = 1,nel
           if (ncell(i,icell) .eq. 2) then
                write(10,*) xcont(i,1+last(i)),
     $                      zcont(i,1+last(i))
                write(10,*) xcont(i,2+last(i)),
     $                      zcont(i,2+last(i))
                write(10,*)
                last(i) = last(i) + 2
           else if (ncell(i,icell) .eq. 4) then
                write(10,*) xcont(i,1+last(i)),
     $                      zcont(i,1+last(i))
                write(10,*) xcont(i,2+last(i)),
     $                      zcont(i,2+last(i))
                write(10,*) xcont(i,3+last(i)),
     $                      zcont(i,3+last(i))
                write(10,*) xcont(i,4+last(i)),
     $                      zcont(i,4+last(i))
                write(10,*) xcont(i,1+last(i)),
     $                      zcont(i,1+last(i))
                write(10,*) xcont(i,3+last(i)),
     $                      zcont(i,3+last(i))
                write(10,*) xcont(i,2+last(i)),
     $                      zcont(i,2+last(i))
                write(10,*) xcont(i,4+last(i)),
     $                      zcont(i,4+last(i))
                write(10,*)
                last(i) = last(i) + 4
           end if
  650      continue
  700 continue
      return
      end
c
      subroutine wrbdy(x,z,nbot,npts,lbot)
      integer npts, lbot
      real x(npts), z(npts)
      integer nbot(0:lbot)
c     Write the coordinates of the points on the boundary to another
c     file.
      do 900 j=1,nbot(1)-1
           write(12,*) x(j),z(j)
  900 continue
      do 1000 jj=1,lbot
           j=nbot(jj)-1
           write(12,*) x(j),z(j)
 1000 continue
      do 1100 j=nbot(lbot)-1,nbot(lbot-1),-1
           write(12,*) x(j),z(j)
 1100 continue
      do 1200 jj=lbot-1,0,-1
           j=nbot(jj)
           write(12,*) x(j),z(j)
 1200 continue
      return
      end
c
      subroutine wrlabl(xcont,zcont,value,len,nlvl,mxlen)
      integer nlvl, mxlen
      real xcont(nlvl,mxlen), zcont(nlvl,mxlen), value(nlvl)
      integer len(nlvl)
c     Write the labels for the contours to a file for use by gnuplot
      write(13,*) 'set parametric'
      write(13,*) 'set data style lines'
      ii = 10
c     (allow labels 1 through 9 for other information)
c     If necessary, remove some of the labels manually or add new ones
c     Modify the labels manually to remove blank spaces and superfluous
c     zeros.
      do 1400 i=1,nlvl
           if (len(i) .ne. 0) then
c               write the label once, in the middle of the contour
                write(13,1300) ii,value(i), xcont(i,len(i)/2 + 1),
     $                         zcont(i,len(i)/2 + 1)
 1300           format('set label ',i3,' "',f10.3,'" at ',f10.3,
     $                 ',',f10.3,' center')
                ii = ii + 1
           end if
 1400 continue
      return
      end
